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We use quantum Monte Carlo methods to demonstrate that the quantum phase diagram of the S=l Heisenberg 
model with uniaxial anisotropy contains an extended supersolid phase. We also show that this Hamiltonian is 
a particular case of a more general and ubiquitous model that describes the low energy spectrum of a class of 
isotropic diwd frustrated spin systems. This crucial result provides the required guidance for finding experimental 
realizations of a spin supersolid state. 
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Theoretical proposals ^ for studying the Bose-Einstein 
condensation (BEC) with magnetic systems were followed by 
a vast number of experimental works . These studies were 
done on spin dimer compounds for which the relevant U(l) 
symmetry is "protected" by the intra-dimer singlet- triplet spin 
gap Isl] . Magnetic systems have the advantage that the mag- 
netic field, which plays role of the chemical potential, can be 
varied continuously over a large range of values. A natural 
question that arises is whether other phases that have been 
proposed for bosonic gases of atoms can be realized in quan- 
tum magnets. The supersolid (SS) state is a prominent and in- 
teresting example because the experimental evidence for this 
novel phase is still inconclusive [4] . 

The search for the SS phase has motivated the study of 
different models for hard-core bosons on frustrated lattices 
111]. These models are relevant for gases of atoms in a peri- 
odic potential (substrate or optical lattice). However, the spin 
S = 1/2 Hamiltonians that are obtained from these models 
by applying a Matsubara-Matsuda transformation [6] are not 
relevant for real magnetic systems. What makes these mod- 
els unrealistic for magnetic systems is the large uniaxial ex- 
change anisotropy. Moreover, the longitudinal and the trans- 
verse components of the exchange interaction have opposite 
signs: while the Ising interaction is antiferromagnetic (AFM), 
the transverse exchange coupling is ferromagnetic. It is then 
natural and relevant to ask if a SS spin phase can exist in a 
magnetic system with isotropic (Heisenberg) interactions. In 
this letter, we provide an affirmative answer to this question 
by calculating the quantum phase diagram of an 5 = 1 spin- 
dimer Heisenberg model. The spin SS phase is induced by 
the application of a magnetic field whose Zeeman splitting is 
comparable to the magnitude of the exchange interactions. 

To understand the physical origin of the spin SS, we shall 
start by considering the simplest (although non-realistic) S = 
1 Hamiltonian that contains this phase in its phase diagram. 
This is an 6* = 1-Heisenberg model with uniaxial single-ion 
and exchange anisotropics on a square lattice: 

Hh = jY^{sfs^^sysy^/\s?si)^Y^{Dsf-Bs?) 

(iJ> i 

(1) 

where (i, j) indicates that i and j are nearest neighbor sites. 



D is the amplitude of the single ion-anisotropy and A de- 
termines the magnitude of the exchange uniaxial anisotropy. 
Note that although the exchange interaction is anisotropic, the 
longitudinal ( J) and transverse (A) couplings are both AFM 
(positive). Henceforth, J is set to unity and all the parameters 
are expressed in units of J. 

The quantum phase diagrams for the spin models consid- 
ered in this paper were obtained by using the Stochastic Se- 
ries expansion (SSE) quantum Monte Carlo (QMC) method. 
The simulations were carried out on a square lattice of size 
N = L X L, with 8 < L < 16. We find rapid convergence 
with N for the system sizes studied (see Fig. [T]). Since the 
SSE is formulated in the grand canonical ensemble, the simu- 
lations are performed at fixed magnetic field, instead of fixed 
magnetization. 

As the external field, B, is varied, the ground state of Hh 
goes through a number of phases, including spin-gapped (IS) 
phases with Ising-like ordering and gapless (XY) phases with 
dominant XY-ordering. The IS phases are characterized by 
long-range (staggered) diagonal order measured by the longi- 
tudinal component of the static structure factor (SSE), 



5-(q) 



-iq-(rj 



(2) 



The XY-phase has long-range off-diagonal ordering mea- 
sured by the transverse component of the SSF, 



iq-(rj-rfc)^^+ 



J 



(3) 



The XY-ordering is equivalent to a Bose-Einstein condensa- 
tion (BEC) whose condensate fraction is equal to 5'+~(Q) 
with Q being the ordering wave-vector (Q = (tt, tt) for the 
case under consideration). The superfluid density corresponds 
to the spin stiffness, ps, defined as the response of the system 
to a twist in the boundary conditions. The stiffness is obtained 
from the winding numbers of the world lines (Wx and Wy) in 
the X- and y- directions: ps = {W^ + Wy) /2(3. 

The IS (XY) phase is marked by a diverging value of 
S^^(Ct) (X N (5'+-(Q) (X TV) in the thermodynamic limit 
A/" ^ oo. In addition, ps vanishes in the gapped IS phase 
while it is finite in the gapless XY phase. A spin SS phase 
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is characterized by a finite value of both S^^{Q)/N and pg. 
Both quantities are always finite for finite size systems and 
estimates for N ^ oo are obtained from finite- size scaling. 

Fig. [T] shows the quantum phase diagram as a function of 
magnetic field, B, for D = 1.5 and A = 1.8. For clarity of 
presentation, S^^{Q) and ps are plotted as a function of the 
resulting magnetization m^. The mz{B) curve features two 
prominent plateaus corresponding to different IS phases. For 
small B, the ground state is a gapped AFM solid (ISl) with no 
net magnetization. The stiffness, ps, and (Q) vanish in 
the thermodynamic limit, while S^^{Q) /N is slightly smaller 
than 1 because the spins are mainly in the S'f = ±1 states 
depending on which sublattice they belong to. The magne- 
tization stays zero up to the critical field, Bd, that marks a 
second order BEC quantum phase transition (QPT) to a state 
with a finite fraction of spins in the = state. This state 
has a finite S^^{Q)/N as well as finite ps and 6'+-(Q)/7V, 
i.e., SS order. The diagonal order results from the 6'f = ±1 
sublattices while the off-diagonal order arises out of a BEC 
of the flipped spins (S^ = "particles"). The magnetization 
increases continuously up to 5 ^ 6.4, where there is a sec- 
ond BEC-QPT to a second Ising-like state (IS2) where all the 
S? = -1 have been flipped to the S? = state. 5'^^(Q)/7V 
remains divergent for TV ^ oo, but the stiffness, ps, and con- 
densate fraction, 6'+~(Q), drop to zero. The ground state 
remains in the IS 2 phase for 6.4 < B < 7.2. Upon further in- 
creasing the field, there is a first order transition to a pure XY- 
AFM phase ( changes discontinuously from = 0.5 to 
rriz ^ 0.59). In the grand canonical ensemble, no ground 
state with any intermediate value of the magnetization is re- 
alized. For a canonical ensemble with a fixed magnetization 
—0.6 < rriz < —0.5, the ground state will phase separate 
into IS2 and XY regions with rriz = 0.5 and rriz = 0.59. In 
the pure XY phase, the diagonal order vanishes while ps and 
(Q)/A remain finite. This situation persists until all the 
spins have flipped to the 6'f = 1 (fully polarized) state. 

Further insight into the SS phase is obtained from the mo- 
mentum dependence of 6'+~(q) ( Fig. [He)). The peaks at 
q = (0,0) and q = Q indicate that the off-diagonal long 
range order is modulated by the presence of solid order. This 
confirms that the SF component of the SS phase results from 
a BEC of S? = spin states that occupy the = —1 or 
down-sublattice with higher probability. This feature distin- 
guishes the SS phase from a uniform canted AFM phase. In 
experiments with magnetic materials, the components of the 
structure factor are selectively measured using standard polar- 
ized neutron scattering experiments. 

For smaller values of A(< D), the second magnetization 
plateau disappears completely (Figl2| leaving a second order 
transition from the SS to the XY-phase. Consequently, there 
is no phase separation regime. The extent of the SS phase 
decreases with decreasing A and vanishes for A 1. 

We shall now discuss the relevance of these results for find- 
ing a SS phase in real magnets. We note that although a 
U(l) invariant model provides a good description of spin com- 
pounds whose anisotropy terms are very small compared to 




FIG. 1 : (Color online) Quantum phase diagram of Hh (Eq. [B for 
D = 1.5 and A = 1.8. (a) Magnetization as a function of field B. 
The SS phase appears between the two Ising (or solid) orderings de- 
noted by ISl and IS2. At higher fields, there is a first order transition 
between the IS2 and the pure XY-AFM phases, (b) Square of the 
XY-AFM order parameter as a function of B.(c) and (d) Longitu- 
dinal component of the staggered SSF and stiffness as a function of 
the magnetization. In a grand canonical ensemble, no ground state 
with 0.5 < rriz ^0.59 (marked PS) is realized - this corresponds to 
the discontinuous IS2-XY transition. For a canonical ensemble with 
magnetization in this range, the ground state phase separates into 
spatial domains with rriz = 0-5 and rriz ~ 0.59. (e) Full momen- 
tum distribution of the form factor, (q). The peak at q = Q, in 
addition the one at q = indicate that the off-diagonal order is mod- 
ulated by the presence of simultaneous long-range diagonal order. 



the Heisenberg interactions, this invariance is never perfect. 
The transition metal magnetic ions belong to this class be- 
cause the spin-orbit interaction is much smaller than the crys- 
tal field splitting. These spin systems have small exchange 
anisotropics for the same reason. Therefore, models that as- 
sume opposite signs for Jz and J_\_ or large values values 
of J±/Jz Mi are not applicable to these spin compounds. We 
will show below that it is not necessary to assume a strong 
uniaxial exchange anisotropy for obtaining a SS phase. 

The system to be considered is a square lattice of 5 = 1 
dimers (Fig. [3]) which only includes isotropic (Heisenberg) 
AFM interactions, an intra-dimer exchange Jo and inter-dimer 
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FIG. 2: (Color online) Same as Fig. [T] but with parameters D = 
1.5, A = 1.2. The second magnetization plateau disappears com- 
pletely. Instead, there is a direct (continuous) SS-XY transition. At 
high fields, there is a transition to a fully polarized state (PL). 



V = AJ, fi = D ^ B 8ind tj = ^j^^ j = 1,2,3 

after we map on each site the eigenstates of onto the eigen- 
states of nil S? = n-, - I and = gl[V2 ^ (1 - \/2)ni]. 
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FIG. 3: (Color online) Square lattice of S=l dimers with an intra- 
dimer Heisenberg AFM interaction Jo and inter-dimer interactions 
Ji and J2. The left side shows the low energy subspace of the single 
dimer spectrum in the presence of a magnetic field. 



frustrated couplings Ji and J2 : 

Hd = -^0 ^ Si+ • Si_ + Ji ^ Sia • Sja 

+ J2 ^ Si«-Sj^-5^5f,. (4) 

(i,j),Q! ia 

The index a = ± denotes the two spins on each dimer. The 
single dimer spectrum consists of a singlet, a triplet and a 
quintuplet (see Fig. [3]). The energy difference between the 
singlet and the triplet is Jo, while the difference between the 
quintuplet and the triplet is 2 Jq. 

For Ji, J2 <C Jo, the low energy subspace of Hd con- 
sists of the singlet, the = 1 triplet and the = 2 
quintuplet (see Fig. [3]). The low energy effective model, H, 
that results from restricting Hd to this subspace is conve- 
niently expressed in terms of semi-hard-core bosonic oper- 
t 

ators, gl and g^, that satisfy the exclusion condition gl = 
(no more than two per site) |H, [13] and obey the commuta- 
tion relations of canonical bosons except for the commutator 
= (5ij(l — ni) (rii = gjg^ is the number operator). The 
expression of H in terms of these operators is: 

(iJ> i 

+ |^ni(ni-l) + F^(ni-l)(nj-l) (5) 

i (iJ> 

with hi = ti(nij - 2)(nij - 3), /12 = 2t2(nij - 1)(3 - nij), 
hs = ts{nij — l)(nij — 2), and = ni + nj. The ampli- 
tudes ti, t2 and ^3 correspond to single-particle hopping terms 
when there are one, two or three particles respectively on the 
corresponding bond The case ti = t2 = ts = t corre- 
sponds to the bosonic Hubbard model with n.n. repulsion JtI] 
in a truncated Hilbert space. Our S = 1 Heisenberg Hamil- 
tonian with uniaxial anisotropy, Hh, is obtained for /7 = D, 



As we mentioned before, H also describes the low energy 
spectrum of Hd- In this case, we have U = Jo,V = (Ji + 
J2)/2, 11 = B - Jq - z{Ji + J2)/2 and tj = 8a^{Ji - 
^2)73 a/3 with a = a/3/2, after mapping the eigenstates of 
+ into the eigenstates of n{ by the simple relations: 

m = + (see FigE) and gl = (^+ - S^_)[^ + 

(1 — ^^)^f]- Fig. [4] shows the quantum phase diagram as 
a function of fi (or B) for = 30.0 and V = 7.0 (in this 
case we take (Ji — J2)/2 as the unit of energy). This set 
of parameters corresponds to Jq = 30 Ji = 8 and J2 = 6 
that satisfies the conditions Jq > z{Ji -\- J2)/2 and Jq ^ 
z{Ji — J2)/2 necessary for the validity of Jf as a low energy 
effective model for Hd- 

At small fi or B, the empty state (all the dimers in a sin- 
glet state) has the lowest energy. For /i > fici (B > Bd) a 
finite density of bosons (triplets) is stabilized in the ground 
state giving rise to a BEG (XY-AFM ordering) at T = 
with a finite the stiffness ps - The absence of solid (Ising) or- 
dering is indicated by S^^{Q)/N 0. The density (mag- 
netization) increases monotonically as a function of /i or 5 
until /i = ^ 2.9 where there is a discontinuous transi- 
tion to a charge-density wave (CDW) or Ising-like phase with 
n = rriz = O.b (the dimers of one sublattice are in a triplet 
state while the other dimers remain in the singlet state). For 
M > Mc3 ^ 23 A, some of the dimers of the singlet sublattice 
are turned into triplets that propagate primarily on the singlet 
sublattice (U > zV where z = 4 is the coordination number). 
Consequently, there is a BEC-QPT (broken U(l) symmetry 
under rotations around the z-axis) inV = d -\-2 dimensions 
to a SS phase, where d is the spatial dimensionality. The diag- 
onal or solid order disappears at an Ising-like quantum critical 
point inV = d -\- 1 dimensions for /i = /ic4 ^ 25.4 (broken 
Z2 symmetry of translation by one lattice parameter followed 
by a TT-rotation around the z-axis). Upon further increase in 
/i, the filling increases monotonically in the resulting SF phase 
till the ground state enters a Mott insulating (MI) phase with 
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FIG. 4: (Color online) Quantum phase diagram of for U = 
30.0, V = 7.0. (a) Particle density n or ruz as a function of the 
chemical potential /x (lower axis) or field B (upper axis), (b) Con- 
densate fraction or square of the AFM-XY order parameter, (c) and 
(d) The staggered SSF and stiffness as a function of n = m^. The 
range of densities marked PS is inaccessible in the grand canonical 
ensemble and would result in a phase separated state in a canonical 
ensemble. 



all the dimers in the triplet state. 

The mechanism for the formation of the SS phase is ex- 
plained most readily in the bosonic language In the 
strong coupling limit (U^V t), the half-filled ground state 
(n = ^) is a checkerboard solid (one sublattice is single occu- 
pied while the other sublattice is empty). Doping away from 
n = ^ results in different scenarios depending on the nature 
of doping and the relation between the coupling constants U 
and V. Extracting bosons from the n = 1/2 crystal costs 
chemical potential /i but no potential energy. The kinetic en- 
ergy gain of the resulting holes is quadratic in t for isolated 
holes (0{t'^/V)), but becomes linear in t if the holes segre- 
gate in a SF bubble. Consequently, if the total density is fixed, 
the system separates in a commensurate crystal with n = 1/2 
and and a uniform SF region with n < 1/2. This implies a 
first order transition between the solid and the SF phases as a 
function of /i (see Figs. [BandlH). 

Doping of the n = 1/2 crystal with additional bosons 
works differently depending on the relation between V and U. 
The energy cost to place a boson at an empty (occupied) site 
is Eq = zV — fi (El =U — fj). Respectively, for U > zV, 
the additional bosons fill empty sites and mask the checker- 
board modulation; foxU — zV ^ |t| the situation is precisely 
particle-hole conjugate to hole doping. In particular, in the 
hard-core limit U ^ oo, the crystalline order is always un- 
stable for n 7^ 1/2. However, for zV ~ U, the bosons can 



be placed on either an occupied or unoccupied site. The ki- 
netic energy gain of the added boson is now linear in t be- 
cause the potential barrier, \zV — U\,fov moving the bosons 
to nearest neighbors is not much bigger than t. As a result, 
the added bosons form a SF phase on top of the density wave 
background and hence the ground state has simultaneous solid 
and SF orders. This SS phase is stable for a sufficiently small 
concentration of added bosons. This is confirmed by the quan- 
tum phase diagram shown in Fig |4] where the SS phase appears 
right next to the n = 1 /2 CDW. We emphasize that this phase 
requires to have two bosons on the same site, which is not pos- 
sible for hard-core bosons (or, equivalently, for 5* = ^ spins). 

In summary, we have shown that simple two-dimensional 
S = 1 Heisenberg models have a spin SS ground state in- 
duced by magnetic field. The physical mechanism that leads 
to this phase does not depend on the dimensionality and simi- 
lar results are expected for three and one-dimensional lattices 
lllll Il2|] . In particular, by finding a SS phase in an isotropic 



8=1 Heisenberg model we are providing the required guid- 
ance for finding this novel phase in real spin systems. The 
crucial ingredients for the described mechanism are: S = 1 
spins, a dimerized structure and sl frustrated inter-dimer cou- 
pling. 
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